clear;clc;

load Factors yyyymm;
TBeg = find(yyyymm==201112); TPort = find(yyyymm==201106); 

load CRSPMonthly permno Return Price ShrOut Volume Excd N T;
ME = Price.*ShrOut/1000; Turn = Volume./ShrOut; clear Price ShrOut Volume;
Return11 = NaN(N,T,'single');
for t = 11:T
    Return11(:,t) = prod(1+Return(:,t-10:t),2)-1;
end
Micro = false(N,T);
for t = 1:T
    idx = Excd(:,t)==1;
    x = ME(idx,t); x(~isfinite(x)) = []; x = sort(x);
    Ns = length(x)/10;
    BrkPts = x([1 floor(Ns*(1:10))]);
    idx = ME(:,t)<=BrkPts(3);
    Micro(idx,t) = true;
end
clear t idx x Ns BrkPts;

load BookValue BE_Monthly; BE = BE_Monthly; clear BE_Monthly;
BM = BE./ME;
load CCMAnnualMonthly AssetsGrowth OperatingProfits;
Prof = OperatingProfits./BE; clear OperatingProfits;
Invst = AssetsGrowth; clear AssetsGrowth;

load Data_CRSPMonthly CommonSampleMonthly;
CommonSample = false(N,T);
CommonSample(:,TBeg:T) = CommonSampleMonthly;
CommonSample(:,TPort:TBeg) = repmat(CommonSampleMonthly(:,2),1,TBeg-TPort+1);
clear CommonSampleMonthly;

NumPort = 10;
temp = struct('AnomalyX',[],'RetPort',[],'WtPort',[]);
VarNames = {'ME', 'BM', 'Prof', 'Invst', 'Return11', 'Return'};

load Factors yyyymm;
ValidStocks = (Excd==1 | Excd==3) & CommonSample; % use all stocks
Port = repmat(temp,6,1);

for loop = 1:2
    if loop==1
        KeepIndx = [3 9];
    else
        KeepIndx = [6 6];
    end

    for vari = 1:4
        Port(vari).AnomalyX = VarNames(vari);
        AnomalyX = eval(VarNames{vari});
        RetPort = NaN(T,NumPort,2,'single');
        WtPort = zeros(N,T,NumPort,2,'single');
        for t = TPort:12:T
            X = AnomalyX(:,t);

            % change signa of some sorting variables so that hi is always expected high return
            if vari==1 || vari==4
                AnomalyX = -AnomalyX;
            end

            idx = isfinite(X) & Excd(:,t)==1 & isfinite(ME(:,t)) & ValidStocks(:,t); % this is the index of stocks to use in calculating breakpoints
            Xsort = unique(sort(X(idx)));  % use unique to do the correct sorts
            NinPort = length(Xsort)/NumPort;
            BrkPts = Xsort([1 floor(NinPort*(1:NumPort))]);
            BrkPts(1) = -Inf; BrkPts(end) = Inf; % set first breakpoint to -Inf; there can be stocks in portfolios that have X1 lower than the lowest X1 from brekpoint stocks

            % Decile 1
            n = 1;
            for s = t+1:min(t+12,T)
                idxKeep_s = X<BrkPts(KeepIndx(1)) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
                idxNew_s = X<BrkPts(2) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
                % this is the index of stocks to be included in the portfolio in this month
                if s==TPort+1
                    idx = idxNew_s;
                else
                    idx = idxNew_s | (idxOld1_sm1 & idxKeep_s);
                end
                idxOld1_sm1 = idx;
                r = Return(idx,s); me = ME(idx,s-1);
                WtPort(idx,s-1,n,1) = 1/sum(idx); RetPort(s,n,1) = mean(r);
                WtPort(idx,s-1,n,2) = me/sum(me); RetPort(s,n,2) = sum(r.*me)/sum(me);
            end

            % Decile 10
            n = 10;
            for s = t+1:min(t+12,T)
                idxKeep_s = X>=BrkPts(KeepIndx(2)) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
                idxNew_s = X>=BrkPts(10) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
                % this is the index of stocks to be included in the portfolio in this month
                if s==TPort+1
                    idx = idxNew_s;
                else
                    idx = idxNew_s | (idxOld10_sm1 & idxKeep_s);
                end
                idxOld10_sm1 = idx;
                r = Return(idx,s); me = ME(idx,s-1);
                WtPort(idx,s-1,n,1) = 1/sum(idx); RetPort(s,n,1) = mean(r);
                WtPort(idx,s-1,n,2) = me/sum(me); RetPort(s,n,2) = sum(r.*me)/sum(me);
            end
        end
        Port(vari).RetPort = RetPort(TBeg:end,:,:);
        Port(vari).WtPort = WtPort(:,TBeg:end,:,:);
    end

    vari = 5;
    Port(vari).AnomalyX = VarNames(vari);
    AnomalyX = eval(VarNames{vari});
    RetPort = NaN(T,NumPort,2,'single');
    WtPort = zeros(N,T,NumPort,2,'single');
    for t = TPort:T-1
        X = AnomalyX(:,t-1); % One extra lag
        s = t+1;

        idx = isfinite(X) & Excd(:,t)==1 & isfinite(ME(:,t)) & ValidStocks(:,t); % this is the index of stocks to use in calculating breakpoints
        Xsort = unique(sort(X(idx)));  % use unique to do the correct sorts
        NinPort = length(Xsort)/NumPort;
        BrkPts = Xsort([1 floor(NinPort*(1:NumPort))]);
        BrkPts(1) = -Inf; BrkPts(end) = Inf; % set first breakpoint to -Inf; there can be stocks in portfolios that have X1 lower than the lowest X1 from brekpoint stocks

        % Decile 1
        n = 1;
        idxKeep_s = X<BrkPts(KeepIndx(1)) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
        idxNew_s = X<BrkPts(2) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
        % this is the index of stocks to be included in the portfolio in this month
        if s==TPort+1
            idx = idxNew_s;
        else
            idx = idxNew_s | (idxOld1_sm1 & idxKeep_s);
        end
        idxOld1_sm1 = idx;
        r = Return(idx,s); me = ME(idx,s-1);
        WtPort(idx,s-1,n,1) = 1/sum(idx); RetPort(s,n,1) = mean(r);
        WtPort(idx,s-1,n,2) = me/sum(me); RetPort(s,n,2) = sum(r.*me)/sum(me);

        % Decile 10
        n = 10;
        idxKeep_s = X>=BrkPts(KeepIndx(2)) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
        idxNew_s = X>=BrkPts(10) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
        % this is the index of stocks to be included in the portfolio in this month
        if s==TPort+1
            idx = idxNew_s;
        else
            idx = idxNew_s | (idxOld10_sm1 & idxKeep_s);
        end
        idxOld10_sm1 = idx;
        r = Return(idx,s); me = ME(idx,s-1);
        WtPort(idx,s-1,n,1) = 1/sum(idx); RetPort(s,n,1) = mean(r);
        WtPort(idx,s-1,n,2) = me/sum(me); RetPort(s,n,2) = sum(r.*me)/sum(me);
    end
    Port(vari).RetPort = RetPort(TBeg:end,:,:);
    Port(vari).WtPort = WtPort(:,TBeg:end,:,:);

    vari = 6;
    Port(vari).AnomalyX = VarNames(vari);
    AnomalyX = eval(VarNames{vari});
    AnomalyX = -AnomalyX; % change signs of REVERSAL so that hi is always expected high return
    RetPort = NaN(T,NumPort,2,'single');
    WtPort = zeros(N,T,NumPort,2,'single');
    for t = TPort:T-1
        X = AnomalyX(:,t);
        s = t+1;

        idx = isfinite(X) & Excd(:,t)==1 & isfinite(ME(:,t)) & ValidStocks(:,t); % this is the index of stocks to use in calculating breakpoints
        Xsort = unique(sort(X(idx)));  % use unique to do the correct sorts
        NinPort = length(Xsort)/NumPort;
        BrkPts = Xsort([1 floor(NinPort*(1:NumPort))]);
        BrkPts(1) = -Inf; BrkPts(end) = Inf; % set first breakpoint to -Inf; there can be stocks in portfolios that have X1 lower than the lowest X1 from brekpoint stocks

        % Decile 1
        n = 1;
        idxKeep_s = X<BrkPts(KeepIndx(1)) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
        idxNew_s = X<BrkPts(2) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
        % this is the index of stocks to be included in the portfolio in this month
        if s==TPort+1
            idx = idxNew_s;
        else
            idx = idxNew_s | (idxOld1_sm1 & idxKeep_s);
        end
        idxOld1_sm1 = idx;
        r = Return(idx,s); me = ME(idx,s-1);
        WtPort(idx,s-1,n,1) = 1/sum(idx); RetPort(s,n,1) = mean(r);
        WtPort(idx,s-1,n,2) = me/sum(me); RetPort(s,n,2) = sum(r.*me)/sum(me);

        % Decile 10
        n = 10;
        idxKeep_s = X>=BrkPts(KeepIndx(2)) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
        idxNew_s = X>=BrkPts(10) & isfinite(Return(:,s)) & isfinite(ME(:,s-1)) & ValidStocks(:,t);
        % this is the index of stocks to be included in the portfolio in this month
        if s==TPort+1
            idx = idxNew_s;
        else
            idx = idxNew_s | (idxOld10_sm1 & idxKeep_s);
        end
        idxOld10_sm1 = idx;
        r = Return(idx,s); me = ME(idx,s-1);
        WtPort(idx,s-1,n,1) = 1/sum(idx); RetPort(s,n,1) = mean(r);
        WtPort(idx,s-1,n,2) = me/sum(me); RetPort(s,n,2) = sum(r.*me)/sum(me);

    end
    Port(vari).RetPort = RetPort(TBeg:end,:,:);
    Port(vari).WtPort = WtPort(:,TBeg:end,:,:);

    yyyymm = yyyymm(TBeg:end);

    if loop==1
        save Portfolios_s20 Port permno yyyymm NumPort -v7.3;
    else
        save Portfolios_s50 Port permno yyyymm NumPort -v7.3;
    end
end
